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Abstract 

This paper detailedly discusses the locally one-dimensional numerical methods for ef- 
ficiently solving the three-dimensional fractional partial differential equations, including 
fractional advection diffusion equation and Riesz fractional diffusion equation. The second 
order finite difference scheme is used to discretize the space fractional derivative and the 
Crank-Nicolson procedure to the time derivative. We theoretically prove and numerically 
verify that the presented numerical methods are unconditionally stable and second order 
convergent in both space and time directions. In particular, for the Riesz fractional dif- 
fusion equation, the idea of reducing the splitting error is used to further improve the 
algorithm, and the unconditional stability and convergency are also strictly proved and 
numerically verified for the improved scheme. 

Mathematics subject classification: 26A33, 65L20. 

Key words: Fractional partial differential equations, Numerical stability, Locally one di- 
mensional method, Crank-Nicolson procedure, Alternating direction implicit method. 



1. Introduction 

The history of fractional calculus can goes back to more than three hundred years ago [12] , 
almost the same as classical calculus. Nowadays it has become more and more popular among 
various scientific fields, covering anomalous diffusion, materials and mechanical, signal process- 
ing and systems identification, control and robotics, rheology, fluid flow, signal processing, and 
electrical networks et al [T3]- Meanwhile, the diverse fractional partial differential equations 
(fractional PDEs), as models, appear naturally in the corresponding field. 

There are already some important progress for numerically solving the fractional PDEs. 
The methods used for classical PDEs are well extended to fractional PDEs, for example, the 
finite difference method [21 EH [SI [2D] > finite element method [H [8] , and spectral method [14] . 
However, almost all of them concentrate on one or two dimensional problems. There are 
already some good developments for realizing the operator splitting (locally one dimension) to 
solve the classical PDEs. This paper focuses on extending the alternating direction implicit 
(ADI) methods to the three-dimensional fractional PDEs, and improving their efficiency. 

The Peaceman and Rachford alternating direction implicit method (PR-ADI) [IB] works 
well for two-dimensional problems. But it can not be extended to higher dimensional problems. 
Douglas type alternating direction implicit methods (D-ADI) [7] are valid for any dimen- 
sional equations. And PR-ADI and D-ADI are equivalent in two dimensional problems. In this 
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paper, we consider the following three-dimensional fractional advection diffusion equation, 

du ^ z ^ =d- XL D~u(x, y, z, t) + d% x D° R u(x, y, z, t) 

+ d i y L D y u ( x , z -> *) + d \ y D^ R u(x, y, z, t) 

+ d i z L D]u(x, y, z, t) + d| z DJ R u(x, y, z, t) (1.1) 
du(x,y,z,t) du(x,y,z,t) du(x,y,z,t) 
ox ay oz 

+ f{x,y,z,t), 

and the Riesz fractional diffusion equation 

, gf J ~ = d i{x L D"u{x,y,z,t) + x D% R u(x,y, z,t)) 

+dt( yi PyU(x,y,z,t) + y D^ R u(x,y 7 z,t)) 
+ d i ( ZL D Zu(x, y, z, t) + z D2 R u(x, y, z, t)) 
+f(x,y,z,t), 

both with the initial condition 

u(x,y,z,0) = u (x,y,z) for (x,y,z)ett, (1.2) 

and the Dirichlet boundary condition 

u(x,y,z,t) = for (x, y, z, t) e 90 x (0, T], (1.3) 

where f2 = (xl,xr)x (yL,yn) X (zl,zr) C I 3 , < f < T, and the fractional orders 1 < a, /3, 7 < 
2; and /(a;, y, z, t) is a forcing function; and all the coefficients are non-negative constants. The 
fractional derivatives used in (jl.lj) and (1.1') are defined as, for 1 < // < 2, 

^x*) = ^-^s r - ( l4 ) 



r(2 - /it) 5a; 2 

and 



= 1X2^) ^ r « - ^-""(ck- (1.5) 

From the viewpoint of conversation law, the advection term in the advection diffusion equation 
should be first order classical derivative, and the fractional derivative corresponding to the 
diffusion term should be Riemann-Liouville one. 

For the two-dimensional case of (|l.ll) - (|1.3p . PR-ADI and D-ADI are discussed and we show 
that they are equivalent for two-dimensional equations. We use D-ADI for the three-dimensional 
(ll.l[) - (|1.3p . The second order finite difference scheme is used to discretize the space fractional 
derivative and the Crank-Nicolson procedure to the time direction. We theoretically prove and 
numerically confirm that the given numerical schemes are unconditionally stable and second 
order convergent in both space and time directions. In general, the ADI methods introduce 
new error term, called the splitting error, comparing with the original discretizations. Usually 
the splitting error term does not affect the convergent order, but most of the time it lowers 
the accuracy seriously. For (1.1'), we use the idea in [7] to reduce the splitting error from 
0(t 2 ) to C(t 3 ) at reasonable computational cost and then recover the accuracy of the original 
discretization, the improved ADI will be called D-ADI-II. The fractional step (FS) method is 
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also simply discussed to show that, after a minor modification to reduce the splitting error from 
0(t) to 0(t 3 ), it is equivalent to D-ADI-II. 

The outline of this paper is as follows. In Section 2, we introduce the second order finite 
difference schemes for the left and right Riemann-Liouville fractional derivatives (|1.4|) and fj 1 . 5|) . 
and the full discretization schemes of the one-dimensional and two-dimensional case of II) - 
(|1.3|) and (|l.l[) - (|1.3p itself are detailedly provided. Section 3 discusses improving the accuracy 
and efficiency of ADI, presents D-ADI-II for (1.1'), and shows that, after a minor modification 
of the FS method, it is equivalent to D-ADI-II. We do the convergence and stability analysis 
for the schemes used in this paper in Section 4. The numerical results are given in Section 5 
and we conclude this paper with some discussions in the last section. 



We use fourth subsections to derive the full discretization of (|l.lj) . and the corresponding 
schemes of (1.1') can be obtained by letting df = d%, d\ = d\, d\ — d%, and k x — n y = k z = 0. 
The first subsection introduces the second order finite difference schemes for the left and right 
Riemann-Liouville fractional derivatives (jl.4[) and (|1.5[) in a finite interval given in [TJ based 
on the idea of [TS]. The second to fourth subsection present the D-ADI schemes for the one- 
dimensional and two-dimensional case of (|l.ll) - (|1.3|) and (|1.1[) - (|1.3I) itself, respectively. 

2.1. Discretizations for the left and right Riemann-Liouville fractional derivatives 

Let the mesh points Xi — x^+iAx, < i < N x , yj = yL+jAy, < j < N y , z m = zl + itiAz, 
< m < N z and t n = nr, < n < N t , where Ax = (x B - x L )/N x , Ay = (y R - y L )/N y , 
Az = (zr — zl)/N z , t = T/Nt, i.e., Ax, Ay and Az are the uniform space stepsizes in the 
corresponding directions, r the time stepsize. For fi £ (1,2), the left and right Riemann- 
Liouville space fractional derivatives (11.31) and (|1.4j) have the second-order approximation op- 
erators 8'^ x ufj m and 8' ll _ x u^ J ■ m , respectively, given in a finite domain [2H!], where ufj m 
denotes the approximated value of u(xi,yj, z m ,t n ). 

The approximation operator of (|1.4p is defined by [TJ [TJ] 



2. Discretization Schemes 




(2.1) 



and there exists 



D X U(X ) yj , z m , t n ) 



x—x 



8' l _ l , +x u{x l ,y j ,z m ,t n ) + O(Ax) 2 , 



(2.2) 



where 



1 



8 fi, + x u { x ij Vj i z rrn tn) 



T{A-^){AxY 




and 




/ = 0, 
1 = 1, 
1 = 2, 



g? = I 6-2 5 -^ + 3 3 ^, 

(I + l) 3 -^ - 4l 3 -^ + 6(1 - lf-f 1 



3 — (.i 



l > 3. 



(2.3) 
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Analogously, the approximation operator of (|1.5[) is described as pQ 
where is defined by (|2.3[) . and it holds that 

t„) + 0(Aa;) 2 , (2.5) 

where 

3u,_x u ( x i> Vji z rrn in) = _ ^t)(Ax)^ ^ < 9l u ( x i+l~liV j > z rmtn) ■ 

In the following, we introduce and list some discrete operators which work for the functions 
of three variables x, y, and z: 

n' ,,™ = u i+l,J,m " t -l,j,m , D „ „ _ D / „ . 

a,x i,j,m 2Ax ' a,£ i,j,m a,x i,j.rm 

-, i+1 

°a, + x u i,j : m— _ a )(Aa;) ct Z^^l u i-l+l,j,mi °a, + x u i,j,m — a l°a, + x u i,j,mi (2.6) 

c/ n \ * a n . r'/ n jx H „.n 

°a,-x^i,j,m rv)( Ax) a ' i+l — °a,-.x^i.j,m a 2^a,-X U i,j.m- 



The discrete operators related to the variable x or y in the above also work for functions of two 
variables x and y, e.g., //,//;;, = ^'f^' 1 ' 3 > S 'a,+xKj = r(4-j)Ax° £mj ff 

Similarly, it is easy to get the one-dimensional and two-dimensioanl case of (|2.ip - (|2.6l) . 



Remark 2.1 QlJ) Denoting U n = [it™, mJ, ■ • • , u^_ 1 ] T , and rewriting \2.1)) and {2.1$ as ma- 
trix forms 5' a +x U n = AU n + b\ and S' a _ x U n = BU n + b 2 , respectively, then there exists 
A = B T . 

2.2. Numerical scheme for ID 

Consider the full discretization scheme to the one-dimensional case of (jl.ll) . namely, 

=d\ XL D%u(x, t) + d% x D« R u{x, t) + k x + /(*> *)■ ( 2 - 7 ) 

In the time direction, we use the Crank-Nicolson scheme. The central difference formula, 
left fractional approximation operator (|2.2[) . and right fractional approximation operator (|2.5I) 
are respectively used to discretize the classical second order space derivative, left Riemann- 
Liouville fractional derivative, and right Riemann-Liouville fractional derivative. Taking the 
uniform time step r and space step Ax, and taking m™ as the approximated value of u(xi,t n ) 
and j™ +1 / 2 = f(xi,t n+ i/2), where £„+i/2 = {t n +in+i)/2, using the one-dimensioanl case of 



(2. 



2T|) - (l2Jl) . we can write (j27Tl) as 

u{xi,t n+1 ) - u(xi,t n ) 1 



diS a:+x u( Xl ,t n+1 ) + dfS a+x u(x 2 ,t n ) + d^8 a _ x u{xi,t n+1 ) 

+ dl5' a ^ x u(xi, t n ) + n x D' a x u(xi, t n+ i) + K x D' a x u(xi, t n ) 
+ f{x l ,t n+1/2 ) + 0(T 2 + (Ax) 2 ), 
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where 



u(x i+ i,t n ) - u(Xj-\,t n ) 

2Ax 



Multiplying (I2.8[) by r, we have the following equation 



\ (d x i8' a +x + d%6' ai _ x + K x D' a x ^ u{xu t n+1 ) 

+ \ (dl8> a+x + d x 2 5' at _ x + k x dQ] u{x u t n ) + rf( Xi ,t n+1/2 ) + R?+\ 



where 



< ct(t 2 + (Aa;) 2 ). 
Therefore, the full discretization of (|2.7I) has the following form 

,n+l 



1 + ~ ( <*i*a. + * + d x 2 5' + k x D')] < + r/; +1/2 . 



(2.9) 



(2.10) 



(2.11) 



We can write (12.111) as 



,n+l 



i+1 



u? + - 
1 2 



r(4-a)(Ax) 
(if 



iV^-i+l 



Z=0 
i+1 



r(4-a)(Ax) 



Ea n+1 , _^ 



21 /■„."+! 



r(4-a)(Ax 



Q 5^ 



d" 2 



1=0 



2Ax 



,n+l\ 



1=0 



r(4 - a) (Ax) 



(=0 



1+1/2 



(2-12) 



For the convenience of implementation, we use the matrix form of the grid functions 

TT n — r-,n „" ,,™ 1 T jpn+1/2 _ r f n+l/2 ,n+l/2 ,n+l/2-|T 

[ 1 ^ 2 ) * * * J N x — lj > — L/l ' / 2 I • ■ • ) JNx-1 1 ' 

therefore, the finite difference scheme (|2.12j) can be rewritten as 



-A 



A 



2 \r(4-a)(Aa;) Q Q T(4-a)(Aa;) a Q 2Aa; 



-B 



jjn+l 



where 



A„ 



I -+ 


If 


(if 


A a - 




2 vr(4 


- a)(Aa;) Q 


f r(4- 




So 














f)o 





... o 


33 


,9 2 v 


9? 


9o 





di 



-At 



-D 



U n +rF n+1/2 , 



(2.13) 



9n x 
9%. 



9n*-2 9n x -3 



9? ffo Q 
92 9? 



and B 



1 
-1 1 

0-10 








... o 



1 

-1 
(2.14) 



G 



2.3. PR-ADI and D-ADI schemes for 2D 

We now examine the full discretization scheme to the two-dimensional case of (jl.lj) . i.e. 



du[X d ^ t] =dl XL D«u(x, y, t) + <F 2 x D« R u(x, y, t) + d\ yL D^u(x, y, t) + d\ y D$ R u(x, y, t) 
du(x,y,t) du(x,y,t) 



~\~ Kx~ 



dx 



dy 



+ f(x,y,t). 



(2.15) 



Analogously we still use the Crank-Nicolson scheme to do the discretization in time direction. 
Taking it? , as the approximated value of u(xi,yj,t n ), using the two-dimensioanl case of (|2.1I) - 
(|2.6|) , we can write (|2.15[) as 



\ \ S $,+y + 6 (S,-y + D 'Lv) u{xi,yj,t n+1 ) 



+ \ (C, + x + <C* + Kx) + \ (Sl +V + 5l_ y + DQ] u(x hyj ,t n ) (2-16) 



+ i-f(xi,y jl t n+1/2 ) + Rlj , 



where 



Kt 1 \<Zr(T J + (Axy + (A V n 



(2.17) 



Using the notations of (|2.6[) . we further dehne 

r .— X 11 l X" 4- n" 
Oa,x ■— ° a , + x T "a,_i "+" U a,. 



Sa.i, ■= S' 



P,V -—°P, + y^°f!,-y + D fS,y 



thus, the resulting discretization of (|2.15D can be written as a Crank-Nicolson type finite dif- 
ference equation 



,,"+1 X n n + 1 _L X <i:™ _L X oi n + 1 _L X a,™ 

U i,p - u i,j _ + Mij + + ,«+l/2 



i.e., 



,n+l 



n+1/2 



(l - - \S^ V ) ujH- 1 = (l + ^ 0)B + + r/, 

The perturbation equation of (|2 . 20[) is of the form 

(i - l§«,x) (i - ~<w) <+* = (i + ~<w) (i + ^,,) «?j + r/;;; l/2 

Comparing (12.211) with (|2.20p , the splitting term is given by 

.n+l „,n 



(2.18) 



(2.19) 



(2.20) 



(2.21) 



(2.22) 



since (w™j~ ~ u ?j) 1S an C( T ) term, it implies that this perturbation contributes an 0[t ) error 



component to the truncation error of the Crank-Nicolson finite difference method (|2.18p . 
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The system of equations denned by (|2.21[) can be solved by the following systems. 
PR-ADI scheme [T5]: 

(i - \& a ,*) < d = (i + \*p,v) <j + ^/;;; l/2 ; (2.23) 

(1 - <r = (1 + <j + \fT /2 - ( 2 - 24 ) 

D-ADI scheme 013 [7]: 

T r- \ ~ ( „ T „ \ „ . „n+l/2. 



1 - J = (1 + ~<W + t<WJ + 17 ; (2.25) 
(l - <f = <j ~ ~<W<r (2-26) 



Take 



Ny-ll U 2,Ny-\T--l U 'N X -l,Ny-\\ ' 

T 



r'i' 



TT™ — I?;" 7/™ 7/™ 7/™ ?; n 

u — L u l,l) u 2,l) • ■ • i a N x — 1,1' "1,2) u 2,2' • ■ • ' a N x — 1,2) ' ' * ' "l, 

r — L/1,1) /2,1) • • • ) JA/a-1,1) /l,2) J2,2) • • • ) JN X -1,2: ■ • ■ > /l,jV v -l) J2,N y -li ■ • ■ J /jV«-l,iV v -lJ ' 

and denote 

g - - 2r(4 - a)(Ai) tt 7 ^ A " + 2r(4 -t) (A*)« f ® ^ + g£ f ® 

^ = 2r(4 - W^VY A "® 1+ 2F(4 - })(Ayy ^ * 1 + 4&j B * *> 

where / denotes the unit matrix and the symbol <g) the Kronecker product |13j . the matrixes 
A a , A/3 and B are defined by (|2.14[) corresponding to a and j3, respectively. Thus, the finite 
difference scheme (I2.21[) has the following form 

(/ - B X )(I - B y )U n+1 = (I + B X )(I + B y )XJ n + tF" +1/2 . (2.28) 

Remark 2.2. The schemes i2.23\) - [K£$ and V2.25)) - \2.26]) are equivalent, since both of them 
come from 112.21)) . 

2.4. D-ADI scheme for 3D 

Using the notations of (|2.1[) - (|2.6I) . we can write (jl.ll) as the following form 

(T T T \ 

1 - ^a,x - -xSp >y - -S ltZ \ u(Xi,yj,Zm,t n+ l) 

( r r r \ +1 ( 2 - 29 ) 

= U + 2$a,x + 2^/3, y + 2 5r V> z ) U ( X ifVh Z m,tn) + Tf{Xi,yj, Z m ,t n+1 / 2 ) + RiJ, m , 

where 



(2.30) 



and 

\Rltl\ ^ ct(t 2 + (Aa;) 2 + (Ay) 2 + (Az) 2 ). (2.31) 





r» 

■ ^aL, + x 


+ S a^x ' 


f D" ■ 




X" 






5~t,z 


■= K+Z 




-D" 
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Similarly, the full discretization scheme of can be written as 



T „ T „ T 



+ 1/2 



(2.32) 



The perturbation equation of (|2.32l) is of the form 



1 - 2^ ) [1- 2 S P-v ) y 1 - 2 5 i' z J u li,m 



(2.33) 



r. 



= (l + -£<*,*) (l + ~<W) (l + ~&r,.) <;, m + fit, 1 ' 2 
The scheme (|2 . 33[) differs from (]2 .32[) by the perturbation term 
(At) 2 (t) 3 

^^(S a , x S 0ty + 5 a . x 5^ z + Sp, v S 7iZ )(u^ m - ul J m ) - ^-5 a , x 8p ty 6 7tZ (<+^ - u^ m ) . 
The system of equations dehned by (12 .33[) can be solved by the D- ADI scheme [5J [7] : 

(l - T -5 a , x ) u$ m = (l + T -5 a , x + tS p , v + rS 7iZ ) u^ m + TfUll 2 ; (2.34) 

(l - \h,v) U i,f,m = U i,i,m ~ \ 5 fS,V u lj,m-> ( 2 ' 35 ) 

~ 2^ 7 ' z ) U h3, m ~ U i,j-rn ~ ~2^1, zU i,j,m- (2.36) 

Similarly, we suppose 

U — L U 1,1,1' U 2,l,l7 • • • ) U N X -X,l,l) ^1,2,1) "2,2, lJ ■ • ■ i u JV a -l,2,li • • • ) 

u l,iVy-l,l> u 2,iV B -l,l) • ' • > u N x -l,N y -l,li 

n n n n n n 

"1,1,2) "2,1,2) ■ • • ) "jVa-1,1,2) "1,2,2) "2,2,2) • • • > "^-1,2,2) • • • ) 

U 1,JV B -1,2) u 2,iV B -l,2) • ■ • ) U N X -1,N V -1,2) ~ * 

■ ■ ■ u l,l,jV«— 1) u 2,l,N z -li ' • ■ ) a N x —l,l,N,—li u l,2,Af z -l) "2,2,iV 2 -l) • • • ) U AT X -1,2,JV Z -1) • • • > ^ 

1 T 

n n n 

"1,AT„-1,JV.-1)"2,JV„-1,JV.-1> • ■ • ) "i\L.-l,Ar„-l,iV.-l ) 



and 



^= 2r(4-T)(A^ / ^ + ^ 

A = WaTT^ 7 ® A P ® 7 + or^ "o'wa...^ 7 A ? ® 7 + TT7. 1 ® B ® 7 > ( 2 - 37 ) 

A = "m/A ~ A ® 7 <g> 7 + ~*' WA — ^ ® 7 <g> 7 + -p^—B <g> 7 ® J, 





2r(4 


- a)(Ax) Q 






2r(4 








2r(4 





2r(4 


- a)(Air) Q 






2r(4 


-/?)(Ay)£ 




^. . A 



AAx 

K v T 
4 Ay' 

4Az' 



where 7 denotes the unit matrix and the symbol <g) the Kronecker product [T3], the matrixes 
Aq,, Ag, A 7 and 7? are defined by (|2.14p corresponding to a, ft and 7, respectively. Thus, the 
finite difference scheme (I2.33[) has the following form 

(I - A X )(I - A y )(I - A z )\3 n+1 = (I + A X )(I + A y )(I + A)U" + rF" +1 / 2 . (2.38) 

The corresponding procedure is executed as follows: 
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(1) First for every fixed z = z m (to = 1, . . . , N z — 1), and each fixed y = yj (j = 1, . . . , N y — 1), 

solving a set of N x — 1 equations denned by (|2.34p at the mesh points Xk,k = 1, . . . ,N X — 1 , 

to get u" j m ; 

(2) Next alternating the spatial direction, and for each fixed x = a;, (i = 1, . . . ,N X — 1), and 

each fixed z = z m (m = 1, . . . , N z — 1), solving a set of iV^ — f equations defined by (|2.35p 
at the points yk, k = 1, . . . , — 1, to obtain u"^ m . 

(3) At last alternating the spatial direction again, and for each fixed y = yj (j = 1, . . . , N y — 1), 

and each fixed x = a;, (i = 1, . . . , iV^ — 1), solving a set of — 1 equations defined by 
(12.361) at the points Zk, k = 1, . . . , N z — 1, to gain w"^- 



3. Improved Accuracy for D-ADI and FS Procedures 

This section shows that the idea of improving the accuracy of D-ADI and FS procedures [7] 
also works well when used to solve Riesz fractional diffusion equation (1.1'). For the simpleness 
to illustrate and discuss this, we focuses on two-dimensional case of (1.1')- It is natural to 
extend higher dimensions. The reason why we abruptly discuss FS procedure here is because 
we want to show FS method is equivalent to D-ADI after some minor modifications even when 
solving fractional PDEs. 



3.1. Correction term for the D-ADI method 

The D-ADI scheme of (12.151) introduces the splitting error term (|2.22p . Even though it is 
still with the order 0(r 2 ), sometimes it will seriously impair the accuracy, see Table l3Tl If we 
add 

to the right hand side of (I2.25[) . then the new D-ADI, called D-ADI-II, will have the splitting 
error 

then the splitting error is reduced to 0(t 3 ). The D-ADI-II is obviously two-step method; ujj 
can be obtained by D-ADI first, then initiate D-ADI-II. 



3.2. Correction term for the FS method 

The original FS method for (|2.15[) should be 



it can be written as 



71,1 _ 

(3.3) 



T 2 a ' xy ^ 

<r - i 



- 2 d f>,y\. u i,j + u i,j)> 



(1 - ^a,«)«# = (1 + 7^K,+rf + 5 ; 

r r ( 3 - 4 ) 
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Comparing (|3.3[) with (|2.20[) , the splitting term is given by 



r 



2 



4 



<W<W<i +u?,), (3.5) 



which is of the order O(r) error component to the truncation error of the Crank-Nicolson finite 
difference method (|2.18p . However, if we add 

T ^5«,^,y{Z<j ~ <,J X ) (3-6) 

to the right hand side of the first equation of p.4[) . then we get the new FS method, called 
FS-II, with the splitting error (|3.2|) . i.e., the splitting error is reduced to C(r 3 ) . 

The FS-II is equivalent to D-ADI-II, since both of them come from the following perturbation 
equation 

' u n+1 



(l - fa, X ) (l - fa,y) 

= (l + fa, X ) (l + fa,y) «?j + T -5 a , x 8p, y (ul 3 - U^ 1 ) + Tfl 



1/2 



i.e., 



U n+1 — ll 71 1 T 



3.3. Accuracy and efficiency of the D-ADI, D-ADI-II, and FS-II methods 

To check the accuracy and efficiency of the D-ADI, D-ADI-II, and FS-II schemes, we consider 
the two-dimensional case of the Riesz fractional equation (1.1'), on a finite domain < x < 
1, < y < 1, < t < 1, with the coefficients df = d\ = 1, and the initial condition u(x, y, 0) = 
sin((2x) A )sin((2 — 2x) 4 )sin((2y) 2 )sin((2 — 2y) 2 ) and the Dirichlet boundary conditions on the 
rectangle in the form u(0,y,t) = u(x,0,t) = and u(l,y,t) = u(x 7 l,t) = for all t > 0. The 
exact solution to this two-dimensional Riesz fractioanl diffusion equation is 

u(x,y,t) = e -< sin{(2x) 4 )sin((2 - 2x) 4 )sin((2y) 4 )sin((2 - 2y) 4 ). 

By the algorithm given in [3] and above conditions, it is easy to obtain the forcing function 
f(x,y,t) at anywhere of the considered rectangle domain with any desired accuracy. 



Table 3.1: The performance of the D-ADI, D-ADI-II, and FS-II methods with Az = Ay = 1/100, and 
the maximum errors (|5.1|l . 



a = 1.9, P = 1.9 


t = 10Ax 


t — 5Ax 


r = 5Aa;/2 


r = Ax 


D-ADI 


3.3496e-02 


7.6895e-03 


2.0624e-03 


6.0826e-04 


FS-II 


2.2638e-03 


1.7249e-04 


2.7765e-04 


3.3166e-04 


D-ADI-II 


2.2638e-03 


1.7249e-04 


2.7765e-04 


3.3166e-04 



From Table [5TT1 we further verify that the D-ADI-II is equivalent to the FS-II method, and 
they may reduce the perturbation error of D-ADI procedure and improve the accuracy. 
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4. Convergence and Stability Analysis 

In the following, we denote by H the symmetric (or hermitian) part of A if A is real (or 
complex) matrix, and || ■ || the matrix 2-norm. 



Lemma 4.1. fJl [TSf The coefficients gf, [i S (1, 2) defined in 
erties 



satisfy the following prop- 



(1) g% = 1,.< = -4 + 2 3 -" < 0, <& = 6 - 2 5 ~" + 3 3 ^; 

(2) 1 > 5 £ > <?£ > 9l > ■ ■ • > 0; 

oo m 

(3) £#=0, £flf<0,m>2. 



Lemma 4.2. Ji7[ p. 28] A real matrix A of order n is positive definite if and only if its sym- 
metric part H = ■ 
of H are positive. 



metric part H = A+ 2 A is positive definite; H is positive definite if and only if the eigenvalues 



Lemma 4.3. ,'17, p. 184] If A e C nxn , let H = be the hermitian part of A, then for 

any eigenvalue A of A, the real part %i(\(A)) satisfies 

■max 

where X m i n (H) and \max{H) are the minimum and maximum of the eigenvalues of H , respec- 
tively. 

Theorem 4.1. Let matrix A a be defined by {2.14% where a £ (1,2), then for any eigen- 
value A of A a , the real part 5ft(A(A Q )) < 0, and the matrix A a is negative definite. Moreover, 
$l(\(diA a + d 2 Al)) < 0, where di,d 2 > 0, d\ + d\ f 0. 



Proof Let H 



A a +A* 



, by (|2.14j) we know 



H = {hi 



2g? 98+9% 93 
98+98 2 5 f gg+gZ g% 



9s 



98+98 2gf 



JS+92 



9n x -i 

9Nx-i 9n x ~2 



9n x -2 9Nx-i 
9n x -2 



93 

2ff? 98+98 
98 + 92 2 3 ? 



(4.1) 



From Lemma I4TT1 it is easy to check that g$ + g% > 0, and the sum of the absolute value of 
the off-diagonal entries on the row i of matrix H is given by 



N x -1 

n= ^2\hi,j\<-9i- 

3—1 ,3¥=i 



According to the Greschgorin theorem [TTJ p. 135], the eigenvalues of the matrix H are in 
the disks centered at hi i, with radius r^, i.e., the eigenvalues A of the matrix H satisfies 



A - h i;i \ = |A - g"\ < n, 
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it implies that X(H) < 0. From Lemma B~2l and |4~3I we obtain that 3?(A(A a )) < and A is 
negative definite. Taking 

H = (*^+^) + (*^+^) T = ( dl + d 2 ) A « + 2 A « = ( dl + d 2 )H, 

similarly, we can prove 5ft(A(diA Q + d^A^)) < 0. 

In the following, we list some properties of the Kronecker Product. 
Lemma 4.4. [13, p. 140] Let A E K mxn , B E R rxs , C E « nxp , and D E W xt . Then 

[A®B){C®D)=AC®BD (eM mrxpi ). 
Lemma 4.5. [13, p. 140] For all A and B, (A ® B) T = A T ® B T . 

Lemma 4.6. [13, p. 141] Let A E R nxn have eigenvalues {A 4 }™ =1 and B E W nxm have eigen- 
values {/ij-J-JLj. Then the mn eigenvalues of A®B are 

\\Hxi • • • > Ai/i m , A2/X1, . . • , A2jU m , . . . , A n /ii . . . , A n /x m . 

Theorem 4.2. Let Ar, .4^ and .A z 6e defined by j2.37\) . Then 

||(/-A)- 1 (/ + ^)||<1, 

where v = x,y, z. 

Proof. From Lemma T4. 5 1 and (|2.37p . we obtain 

A x + Al _ (d? + df)r _ i ^ i ^(A a + Al 



(4.2) 



2 2r(4-a)(Aa;) a 
A y +A T y ^ {d\ + dpT ( A p +A T p 
2 2T(4-P)(Ay)P ® [ 2 



2 2r(4- 7 )(Az)7 ^ 2 j^ 1 ® 1 - 

According to Theorem 14.11 and Lemma l4~6l we know that A v +A^ are negative definite and 
symmetric matrices, where v = x, y, z. Then for any v = (ui, «2j • ■ • , v n ) T E K n , we have 

u T w < w T (7-^)(J-^)w. 

Substituting u and v T by (7 — A v )~ 1 v and w T (7 — „4^) _1 , respectively, leads to 

v T {I-J^)- 1 {I-A u )- 1 v<v T v. 

Then, there exists 



(J- A) =supW = < 1. 

v^o V w J v 

Similarly, we have 

v T {L + Al){I + A v )v < v T (L - Al){I - A v )v. 
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Taking v by (I — A u ) v, then the above equation can be rewritten as 

v T {I - J£)-\I + Al)(I + A V ){I - A^v < v T v. 
From Lemma T4. 41 it is to check that A x , A y and A z commute, then it yields that 
Wil-A^il + A^W = \\(I + A V )(I - A,)- 1 ]] 



,„ , iPV - Al)-^I + Al){I + A„)(I A v )-^v 
sup y T i !• 



4.1. Stability and Convergence for ID 

Theorem 4.3. The difference scheme h2.11\) with a G (1,2) is unconditionally stable. 

Proof. Let u™ (i = 1,2,..., N x — 1; n = 0, 1, . . . , iV t ) be the approximate solution of u™, 
which is the exact solution of the difference scheme (|2.1ip . Putting e" — uf — uf, and denoting 
E n = [e™, 63 , . . • , e% then from (12 . 1 1[) we obtain the following perturbation equation 

{I-M)E n+1 = {I + M)E n , 

where 

T ( d^ d~^ Kr \ 

M = 2 {r(4-a)(Ax)« A + T{A-a){Axr AT + 2Kx~ B ) ' (4 ' 3) 
Denoting A as an eigenvalue of the matrix M, and using (|4.3[) . there exists 

M + M T r(df + df) 

2 ~~ 2r(4-a)(Ax) Q ' 

where is dehned by (|4.1j) and negative definite by the proof Theorem 14. 11 then from Lemma 
1431 we get 5R(A(M)) < 0. 

Note that A is an eigenvalue of the matrix M if and only if 1 — A is an eigenvalue of the 
matrix I — M, if and only if (1 — A) _1 (l + A) is an eigenvalue of the matrix (I — M) _1 (7 + M). 
Since SR(A(M)) < 0, it implies that |(1 - A) _1 (l + A)| < 1. Thus, the spectral radius of the 
matrix (I — Af) _1 (7 + M) is less than 1, hence the scheme (|2.11[) is unconditionally stable. 

Theorem 4.4. Let u(xi, t n ) be the exact solution of \2. 7\ l with a £ (1, 2), and it" be the solution 
of the finite difference scheme i2.11\) . then there is a positive constant C such that 

\\u(xi,t n ) - u?\\ <C(t 2 + (Ax) 2 ), i = l,2,...,JV x -l;n = 0,l,...,JV t . 

Proof. Denoting e" = u(xi,t n ) — uf, and e n = [e™, , • . • , e-y _J T . Subtracting (I2.9p from 
(|2.11[) and using e° = 0, we obtain 

(I - M)e n+1 = (I + M)e" + R n+1 , 

where M is defined by (|4. 3[) . and R n = [R™ , RV^ , ■ • ■ , Kjj _i] T - The above equation can be 
rewritten as 

e n+1 = (I - M)~ l {I + M)e n + (I - M)~ 1 R n+1 , 

and taking the 2-norm on both sides, similar to the proof of the Theorem 14.21 we can show that 
||(7 - M) _1 (J + M)\\ < 1. Then, using \R" +1 \ < ct(t 2 + (Ax) 2 ) in (j2~TUl) . we obtain 

n-1 

||e"| < |(/-M)- 1 (/ + M)e n - 1 || + ||i?"|| < He™- 1 !) + < ^ ||i? fe+1 || < c(r 2 + (Ax) 2 ). 

fc=0 
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4.2. Stability and Convergence for 2D 

Theorem 4.5. The difference scheme i2.21\) with a,j3€ (1>2) is unconditionally stable. 

Proof. Let ufj (i — 1, 2, . . . , N x — 1; j = 1, 2, . . . , N y — 1; n = 0, 1, . . . , Nt) be the approximate 
solution of ufj, which is the exact solution of the difference scheme (|2.21l) . Taking e"^ — 
ufj — ufj , then from (|2.21l) we obtain the following perturbation equation 

(/ - B X )(I - B y )E n+1 = (I + B X )(I + B y )E n , (4.4) 

where B x and B y are given in (I2.27p . and 

xpn r n n n n n n n n n ]T 

rj — L e l,l> e 2,l) • • ■ i e JVa,-l,l> e l,2) £ 2,2j • • • ) e Ar x -l,2J • • - ) l,N y — l ' 2,N y — l ) • • • ) e JV x - 1, JV H - lJ j 

and we can write (|4.4I) as the following form 

E n+i = (J _ _ S x )- X (7 + B x )(7 + B,)E". (4.5) 

Using Lemma 231 it is to check that B x and B y commute, i.e., 



d\r d\r 
2T(4-(3)(A y y ^ + 2L(4-/3)(Ay)^^4A^ J 



B X B„ = = ( ^ + m/A ^ + — ^ Aj, + ~rx~B 



dlT A a + —Jt——A T a + —^—B 



2T(4:-a)(Ax) a 2r(4 - a)(Ax) a a 4Ax 

Then Eq. (|4.5D can be rewritten as 

E" = ((/ - B v )- l (I + B y )) n (((/ - B*)- 1 ^ + B x ))" E° 
From Lemma T4. 5 1 and (I2.27[) . we obtain 

B. + Bj (df + df)r /A a + 4£ 



(4.6) 



2 2L(4-a)(Ax) Q 

B y + B T y = (4+j|)r / A^+Af 
2 2r(4-/3)(Ay)M 2 



<g> J. 



According to Theorem l4.1l and Lemma l4~2l the eigenvalues of Aa \ A " and ^"^^ are all negative 
when a, f3 € (1,2). Let X x and X y be an eigenvalue of matrices B x and By, respectively. 
From Lemma 14.61 we get !R(A X ) < and $t-{\ y ) < 0, then, the eigenvalues of the matrices 
{I-B x )- l {I + B x ) and (/ - By)- 1 (I + By), (1 - A*)"^! + X x ) and (1 - A^-^l + Aj,) are less 
than one. And it follows that ((I - By)- 1 (I + B y )) n and (((/ - B x )~ l (I + B x )) n converge to 
zero matrix [TTl p. 26], as n —> oo. Hence the difference scheme (|2.19| is unconditionally stable. 



Theorem 4.6. Let u(xi, yj,t n ) be the exact solution of \2. 1 5\) with a, /3 G (1,2), and ufj be 
the solution of the finite difference scheme h2.21\) . then there is a positive constant C such that 

\\u(x uyj ,t n ) - ufj] < C(r 2 + (Ax) 2 + (Ay) 2 ), 



where i = 1, 2, . . . , N x - 1; j = 1, 2, . . . , N y - 1; n = 0, 1, . . . , N t . 
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Proof. Taking e"^ = u(xi,yj,t n ) — ufj, and subtracting (|2 . 16[) from (j2.21j) . we obtain 

(/ - B X )(I - B y )e n+1 = (/ + B X )(I + B y )e n + R"+\ (4.7) 
where B x and B y are given in (|2.27[) . and 

n r n n n n n n n n n ] T 

e — L e l,l) e 2,l> • • • > e JV x -l,li e l,2> e 2,2, • ■ • i e AT x -l,2, • • • > e l,JV !/ -l! e 2,AT H -l, • ■ • > e 7V x -l,7V B -lJ ' 

Rn r pn nil pn rm nil pn nil nn nti ] T 

— [-0.1,1,-0.2,1! • • • ) "-N*-!,!) ^1,21 n 2,2i ■ ■ ■ ; ^iVit-l^! • • • ) ^l.iVy-l) n 2,N v -li ■ ■ ■ , it 7V x - l,iV„-lJ 

and < cY(r 2 + (Ax) 2 + (Ay) 2 ) is given in (l2~T71) . 

Since B x and £\, commutes in (|4.6[) . then Eq. (|4.7I) can be rewritten as 

e „+i = (/ _ Bb )-i (/ + B(e )(j _ + By )e« + (/ - - B T ) _1 R" +1 , 

and taking the 2-norm on both sides, similar to the proof of Theorcm l4.21 it can be proven that 
|| {I-B x )-\I + B x ){l-B y )-\I + By)\\ < \\(I-B x )- l (I + B x )\\-\\(I-B y )-\l + B y )\\ < 1 
and -By)- 1 ^ -B^W < ||(/ - ) 1 1 1 ■ ~ B x )- l \\ < 1. Then we get 

n-l 

||e n || < H Rfc+1 H ^ c ( r2 + ( Ax ) 2 + ( A V) 2 )- 

k=0 

By almost the same proof to the theorems of 2D, we can prove the following results for 3D. 

Theorem 4.7. The difference scheme h2.33]) of the fractional convection diffusion equation 
with a,/3,7 G (1,2) is unconditionally stable. 

Theorem 4.8. Let u(xi, yj, z m , t n ) be the exact solution of U.l\) with a, /3,7 G (1,2), and u™j m 
be the solution of the finite difference scheme \2.33\) . then there is a positive constant C such 
that 

\\u{x uyj ,z m ,t n ) - ul jt J\ < C{t 2 + (Ax) 2 + (Ay) 2 + (Az) 2 ), 
where i = 1, 2, . . . , N x - 1; j = 1, 2, . . . , N y - 1; m = 1, 2, . . . , N z - 1; n = 0, 1, . . . , N t . 

4.3. Stability and Convergence of (1.1') by using D-ADI-II and FS-II 

To prove the stability and convergence of D-ADI-II and FS-II for (1.1'), we need the following 
two lemmas. 

Lemma 4.7. p. 396] If P and P + Q are m-by-m symmetric matrices, then 
Xk(P)+X m (Q) <Xk{P + Q)< A fc (P) + Ai(Q), fc = l,2,...,m, 
with eigenvalues Ai(-) > A2(-) > • ■ • > A m (-). 

Lemma 4.8. [10, p. 84] Let the quadratic equation be A 2 — b\ + c = 0, where b and c are both 
real parameters, then all roots satisfy |A| < 1 if and only if \b\ < 1 + c < 2. 

Theorem 4.9. The difference scheme \3. 7| ) corresponding to two-dimensional case of (L.l') 
with a, [3 € (1,2) is unconditionally stable. 



1G 



Proof. For the two-dimensional case of (1.1'), Eq. (|2.27p has the following form 

d x 



g - = 2 r(4-lW ^ ( ^ + ^ ) = ^' where ^ 



B = tL 

Dy 2T(4-0)(Ay)/» 



W^K^¥^ (Aa + ATal 
dy 



{A [3 +Aj)®I= -B v , where B y = —-———{Ap + Aj 3 ) ® I, 

(4.8) 



and and £> y commute, i.e. 



d x 



T(A-a)(Ax) a r(4-j8)(Ay) 

Let u™ ■ (i = 1, 2, . . . , iVa, — 1; j = 1, 2, . . . , JV y — 1; n = 0, 1, . . . , N t ) be the approximate 
solution of it™,-, which is the exact solution of the difference scheme (|3.7[) . Taking e™, = — it™ •, 



then from (|3.7[) we obtain the following perturbation equation 



(J - B B )(J - B y )E n+1 = (7 + B iC )(7 + £ S )E" + S^E" - B^.E"" 1 , 



i.e., 



where 



E 



n+l 



(P + Q)E™ - QE" 



(4.9) 



E 



and 



L e l,l J e 2,l) • • • J e jVa-l,l) e l,2i e 2,2; • ■ • ) e N w -l,2> • • • ) e l,AT„-l) e 2,A r B -li • • • ) ^a,— 1,;V S -1.I 



P = (I - B X )-\I + B X )(I - By)- 1 ^ + B y ), Q = (I- B x )~ l B x {I - By)- x B y . 



Therefore, Eq. (I4.9[) can be rewritten as 



with 



rn+l 



E rl+1 
E" 



yn+1 = M y» 



and M 



P + Q -Q 
I 



From [21, p. 128], we know that the eigenvalues of M are the same as the eigenvalues of L, 
where 

j _ \X k (P + Q) -A fc (Q)' 
I 1 . ' 

then the eigenvalues A of M satisfies 

A 2 - \ k (P + Q)A + A fc (Q) = 0, k = l,2,...,m(m = N x -l-N y -l). 

Similar to the above proof, we know that B x and B y are negative definite and symmetric 
matrices, and the matrix B x B y or B y B x is positive definite and symmetric, it follows that 
Afc(-P -I- Q) and X k {Q) are real numbers, and we have 

|A fc (P)|<l and 0<A fc (Q)<l, 

According to Lemma 14.71 and 14.81 we get 

-1 - A fe (Q) < A fc (Q) + X m (P) < X k (P + Q)< A fe (Q) + A X (P) < X k (Q) + I, 

thus, the difference scheme is unconditionally stable. 
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Theorem 4.10. Let u(xt, yj,t n ) be the exact solution of \2.15}) corresponding to two-dimensional 
case of (\.V ) with a,f3 G (1,2), and ufj be the solution of the finite difference scheme \S. 7| ). 
then there are a positive constant C and some kind of norm \\ \ ■ \ \\ such that 

Mxi^tn) - <J|| < C(r 2 + (Ax) 2 + (Ay) 2 ), 



where i = 1, 2, . . . , N x - 1; j = 1, 2, . . . , N y - 1; n = 0, 1, . 



Proof. For the two-dimensional case ol (1.1'), taking e"j — u(xi,yj,t n ) — u"j, from (|2.16p 
and p.7p . we obtain 



(I - B X )(I - B y )e n+1 = (I + B X )(I + B y )e n + B x B v (e n - e™" 1 ) 



R 



n+l 



(4.10) 



p>" — \p n p n p n p n p n p n p n p n p n ] T 



where B x and B y are given in (|4.8[) . and 

pn r pn pn pn pn pn pn pn pn 

11 — L-"l,l) rt 2,l, • • ■ , "'N,,— 1,1 J -""1,2) "2,2) • • • ' ^JVaj-l^) • • • J "ljJVy-l) "2,^-1, ' 

and li?™^ 1 ! < cr(r 2 + (Ax) 2 + (Ay) 2 ) is given in ([237) . Similarly, taking 

P = (/ - B*)"^! + B X )(J - S^-^/ + B„); 
Q = (J -B,,)" 1 ^; 
S = {I-B x )-\l-B v )-\ 

then, Eq. (I4.10[) can be rewritten as 

v «+i = MV" +NR ,l+1 , 

with 



pn ]3 
' "■Nx-l.Nv—li 



(4.11) 



yn+l _ 



a n+l 



, R 



n+l _ 



R n 



P + Q -Q 
J 



and N = 



S 




such that 



Similarly, we can prove |Afc(N)| < 1, then there exists some kind of norm 
|||M||| < 1, and |||N||| < 1. Taking the norm on both sides of fiTTj) leads to 

n-l 

llle'll < |||V"||| < m k+1 \W < <t 2 + (Ax) 2 + (Ay) 2 ). 

fe=0 



All the theoretical results for the three-dimensional case (1.1') can be obtained by the same 
way of the two-dimensional case of (1.1'). For the briefness of the paper, we omit them here. 



5. Numerical Results 

Here we verify the above theoretical results including convergent order and stability. Intro- 
ducing the vectors U (Ax) = [v,h(xo, t), . . . , Uh(x n , t)] T , where U is the approximated value, and 
u(Ax) — [u(xo,t), . . . ,u(x n ,t)} T , where u is the exact value and the stepsize in space is Ax, 
i.e., Ax = aij+i — Xi, in the following numerical examples the errors are measured by 

\\U(Ax)-u(Ax)\\ 00 , (5.1) 

where || • is the maximum norm. 
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5.1. Numerical results for ID 



Let us consider the one-dimensional two-sided fractional convection diffusion equation (|2.7[) . 
where < x < 1 and < t < 1, with the coefficients df = df = k x = 1, and the forcing function 



x 2 (1 - xf + (Ax 3 - 6x 2 + 2x) + 



r(3) 

r(3 - a) 



(x 2 - a + (1 - xf- a ) 



2 „ F(4) ; (x 3 - a + (1 - xf- a ) + ^ (5) (x 4 -" + (1 - x) 4 - a ) 



r(4-a) v ~ ' Vi ~' 1 ' r(5-a) 

the initial condition u(x, 0) = x 2 (l — x) 2 , the boundary conditions it(0, t) =u(l,t) =0, and the 
exact solution of the equation is u(x, t) = e^ t x 2 (l — x) 2 . 



Table 5.1: The maximum errors (|5.1|l and convergent orders for the scheme H2.11|l of the one-dimensional 
two-side fractional convection diffusion equation (|2.7|l at t=l and At = Ax. 



At, Ax 


a = 1.1 


Rate 


a = 1.5 


Rate 


a = 1.9 


Rate 


1/10 


0.0022 




0.0011 




0.0010 




1/20 


4.5729e-004 


2.2916 


2.6284e-004 


2.0596 


2.5502e-004 


1.9917 


1/40 


1.0712e-004 


2.0939 


6.2954e-005 


2.0618 


6.4257e-005 


1.9887 


1/80 


2.5242e-005 


2.0853 


1.5067e-005 


2.0628 


1.6169e-005 


1.9906 


1/160 


5.9414e-006 


2.0869 


3.6083e-006 


2.0620 


4.0594e-006 


1.9939 



In Table 15.11 we show that the scheme (|2.1ip is second order convergent in both space and 
time. 

5.2. Numerical results for 2D 



Consider the two-dimensional two-sided space fractional convection diffusion equation (I2.15[) , 
on a finite domain 0<x<2,0<y<2,0<i<2, and with the coefficients 



dl = T(3 - a)x a , d%=r(3-a)(2-x) a , k x = -x, 
d\ =T(3 -0)yP, d^=T(3-0)(2-yf, K y = \y, 
and the forcing function 

f(x, y,t)=- Ae- t x 2 y 2 {x - 2)(y - 2)(3xy - 5x - by + 8) 

3 (x 3 + (2 - x) 3 ) 3 (x 4 + (2 - x) 4 ) 



32e-V(2-2/) 2 
32e~ t x 2 (2 -x) 2 



x 2 + {2-xf 
V 2 + (2 - y) 



3 — a 



(3-a)(4-a) 



, t , 3(y 3 + (2-y) 3 ) 3 (y 4 + (2 - y) 4 ) 
3-/3 (3- 0X4-0) 



and the initial condition u(x, y, 0) = 4x 2 (2 — x) 2 y 2 (2 — y) 2 and the Dirichlet boundary conditions 
on the rectangle in the form u(0, y, t) — u(x, 0, t) = and u(2, y, t) = u(x, 2, t) = for all t > 0. 
The exact solution to this two-dimensional two-sided fractional convection diffusion equation is 



u(x,y,t) = Ae- t x 2 (2 - x) 2 y 2 {2 - yf . 
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Table 5.2: The maximum errors (|5.ip and convergent orders for the scheme (I2.25l) - (|2.26|) of the 
two-dimensional two-sided fractional convection diffusion equation (|2.15|) at t — 2 with r = Ax = Ay. 



T, Ax, Ay 


a = 1.1, = 1.2 


Rate 


a = 1.5, P = 1.4 


Rate 


a = l.9,/3 = 1.9 


Rate 


1/10 


0.0133 




0.0116 




0.0109 




1/20 


0.0033 


1.9966 


0.0029 


1.9830 


0.0028 


1.9750 


1/40 


8.3408e-004 


1.9985 


7.4001e-004 


1.9876 


7.0378e-004 


1.9739 


1/80 


2.0877e-004 


1.9982 


1.8612e-004 


1.9913 


1.7900e-004 


1.9752 


1/160 


5.2231e-005 


1.9990 


4.6726e-005 


1.9940 


4.5468e-005 


1.9770 



Comparing Table 15.21 with Table 2 of [T] , we further confirm that the PR- ADI and D- 
ADI arc equivalent for solving two-dimensional equations, since they have the completely same 
maximum error values. Table f5T21 numerically shows that the D-ADI scheme (|2.25|) - (|2.26|) is 
second order convergent and this is in agreement with the order of the truncation error. 

5.3. Numerical results for 3D 

Consider the three-dimensional two-sided fractional convection diffusion equation (jl.ip . on 
a finite domain < x < 2, < y < 2, < z < 2, < t < 2, and with the coefficients 

tff =r(3-a)a; Q , df = T(3 

d? = r(3-/V, 4 = r(3 

rff = r(3- 7 > 7 , d* 2 = T(3 

and the zero Dirichlet boundary conditions on the cube for all t > 0, the exact solution to this 
three-dimensional two-sided fractional convection diffusion equation is 

u(x, y, z, t) = 4e"*x 2 (2 - x) 2 y 2 (2 - y) 2 z 2 (2 - zf . 

According to the above conditions, it is easy to get the forcing function f(x, y, z, t). 



Table 5.3: The maximum errors (|5.1|l and convergent orders for the scheme (|2.34|l - (|2.36|l of the three- 
dimensional two-sided fractional convection diffusion equation (|1.1[) at t = 2 with r = Aa; = Ay = Az. 



T 


a = P = 7 = 1.2 


Rate 


a = 1.4, P = 1.5,7= 1-6 


Rate 


a = P = 7 = 1.9 


Rate 


1/10 


1.2063e-002 




1.3349e-002 




1.5558e-002 




1/20 


3.0047e-003 


2.0053 


3.3242e-003 


2.0057 


3.7859e-003 


2.0390 


1/40 


7.5079e-004 


2.0008 


8.3225e-004 


1.9979 


9.4168e-004 


2.0073 


1/80 


1.8773e-004 


1.9997 


2.0875e-004 


1.9952 


2.3625e-004 


1.9949 



Table [5751 also shows the maximum error, at time t — 2 and r = Ax = Ay = Az, between 
the exact analytical value and the numerical value obtained by applying the D-ADI scheme 
(|2.34[) - (|2.36[) . and the scheme is second order convergent and this is in agreement with the 
order of the truncation error. 



-a)(2-x) a , k x = -x, 
-f3)(2-yf, K v = \y, 
- 7 )(2-z) 7 , k z = ^z, 
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6. Conclusions 

This work provides an algorithm which can efficiently solve three-dimensional space frac- 
tional PDEs. The idea is to solve higher dimensional problem by the strategy of dimension by 
dimension. When realizing the idea, the splitting errors may be introduced, so the techniques 
of diminishing the influences of splitting errors are also discussed. The effectiveness of the 
algorithm is theoretically proved and numerically verified. 
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